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A COMPARISON  OF  THREE  DIFFUSION  MODELS  OF  THE 
UPPER  MIXED  LAYER  OF  THE  OCEAN 


I.  INTRODUCTION 


In  recent  years  there  has  been  increasing  interest  in  modeling  the 


upper  mixed  layer  of  the  ocean.  This  is  the  turbulent  region  near  the 


surface  which  is  strongly  mixed  by  the  wind.  The  turbulence  is  generated 


either  by  wave  action,  by  shear  instability  in  the  wind  driven  surface 


Etaoan  layer,  or  by  convective  instability. 

The  dynamics  of  the  mixed  layer  is  important  because  it  affects 


obtains  most  of  its  energy  from  latent  and  sensible  heat  exchange  with 


the  ocean,  and  these  fluxes  of  heat  are  strongly  dependent  upon  the  air 


sea  temperature  difference.  The  sea  surface  temperature  in  turn  depends 


upon  the  sea  surface  heat  flux,  the  depth  of  the  mixed  layer,  and  the 
rate  at  which  the  mixed  layer  depth  is  changing.  It  is  also  influenced 
by  the  history  of  past  warnings  and  deepenings  which  determine  the  tem- 
perature structure  below  the  mixed  layer  down  to  the  base  of  the  seasonal 
thermocline.  Therefore,  prediction  of  the  storage  of  heat  in  the  upper 


ocean,  the  sea  surface  temperature,  and  the  air- sea  heat  flux  requires 


careful  modeling  of  the  mixing  processes  in  the  upper  layers  cf  the  sea 
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This  work  consists  of  a description  and  comparison  of  the  eddy- 
diffusion  models  of  Munk  and  Anderson  (19^8),  Vager  and  Zilitinkevich 
(19d8),  and  Mel lor  and  Yamada  (Level  2 model,  19 7M . These  models  are 
parameterizaflons  of  the  vertical  fluxes  of  heat  and  momentum  in  a turbu- 
lent, stratified  boundary  layer.  Such  a parameterization  is  required  to 
close  and  complete  the  equations  for  the  conservation  of  heat  and  momentum 
for  a model  of  the  upper  ocean. 

The  three  models  were  chosen  for  comparison  because  they  are 
fairly  simple,  easily  implemented  turbulence  parameterizations . The  Munk- 
Anderson  model  was  derived  in  an  empirical  fashion  from  observations  of 
stratified  turbulence.  The  Vager-Zilitinkevich  and  Mellor-Yamada  models 
are  more  recant  attempts  to  develop  simple  turbulence  models  that  take 
greater  advantage  of  turbulence  closure  theory.  The  models  were  compared 
to  determine  how  greatly  their  predictions  of  the  depth  of  the  mixed  layer 
differ  for  identical  conditions,  and  some  attempt  is  made  to  relate  the 
differences  in  predicted  mixed  layer  depths  to  the  structural  differences 
of  the  models. 

The  eddy  diffusion  models  were  studied  with  a one- dimensional  model 
in  order  to  concentrate  solely  on  the  vertical  processes  occurring  in  the 
upper  ocean.  Therefore,  only  vertical  diffusion  is  included  in  the 
equations  and  horizontal  advection  and  diffusion  are  ignored.  Vertical 
advection  was  also  ignored  on  the  assumption  that  the  vertical  upwelling 
typical  of  the  open  ocean,  which  could  have  been  included  in  the  model, 
ioes  not  play  an  important  role  in  the  dynamics  of  the  mixed  layer  over 
short  time  periods  (days). 
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Two  problems  of  mixed  layer  deepening  were  examined  with  each  of 
Che  three  diffusion  models.  The  first  problem  is  that  of  the  mixed  layer 
depth  achieved  in  an  initially  motionless,  linearly  stratified  upper 
ocean  due  to  the  imposition  of  a constant  wind  stress  and  zero  surface 
heat  flux.  The  second  problem  is  that  of  the  mixed  layer  depth  achieved 
in  an  Initially  motionless,  weakly  stratified  upper  ocean  due  to  the  im- 
position of  a constant  wind  stress  and  a constant,  positive  heat  flux  at 
the  sea  surface.  (A  positive  surface  heat  flux  Implies  a flux  of  heat 
into  the  sea.)  These  two  problems  are  representative  of  two  distinct 
types  of  mixed  layer  deepening  that  might  be  expected  in  Che  real  ocean. 
Denman  (1973)  refers  to  the  first  problem  as  a case  of  wind  dominated 
deepening  of  the  mixed  layer  and  to  the  second  as  a case  of  heat  domi- 
nated deepening. 

Wind  dominated  deepening  is  characterized  by  a small  or  negligible 
sea  surface  heat  flux.  The  deepening  of  the  mixed  layer  is  determined 
by  the  magnitude  of  the  surface  wind  stress  and  the  stability  of  the 
stratification  of  the  upper  ocean.  Turbulence  generated  by  the  wind 
stress  erodes  the  stable  stratification  at  the  base  of  the  mixed  layer, 
and  as  this  water  is  entrained  into  the  mixed  layer,  the  mixed  layer 
deepens . 

Mixed  layer  models  such  as  those  of  Kraus  and  Turner  (1967)  and 
Denman  (1973)  are  based  on  the  premise  that  a certain  fixed  fraction  of 
the  kinetic  energy  being  imparted  to  the  ocean  by  the  surface  wind  stress 
is  used  to  do  work  on  the  upper  ocean  by  increasing  its  potential  energy 
through  turbulent  mixing.  In  the  case  of  wind  dominated  mixing  where 
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the  surface  heat  flux  is  very  small,  the  wind  stress  can  only  do  work  on 
the  upper  ocean  by  deepening  the  mixed  layer.  As  the  mixed  layer  deep- 
ens, the  potential  energy  of  the  upper  ocean  is  increased  by  the  entrain- 
ment and  upward  transport  of  denser  water  by  turbulence.  Hence,  for  the 
case  of  wind  dominated  mixing,  the  Denman  and  Kraus  and  Turner  models 
predict  that  the  mixed  layer  will  continue  to  deepen  as  long  as  the  wind  ' 
is  blowing.  Observations  of  the  actual  behavior  of  the  ocean's  mixed 
layer  seem  to  indicate  that  this  behavior  is  unrealistic. 

The  models  to  be  compared  here  make  no  assumptions  regarding  the 
fraction  of  kinetic  energy  input  by  the  wind  that  acts  to  increase  the 
potential  energy  of  the  upper  ocean.  It  will  be  seen  that  the  Vager- 
Zllitinkevich  and  Mellor-Yamada  models  predict  that,  for  wind  dominated 
deepening  due  to  the  imposition  of  a constant  wind  stress,  the  deepening 
of  the  mixed  layer  virtually  ceases  -within  a couple  of  days.  At  this 
point,  all  the  kinetic  energy  generated  by  the  wind  is  lost  to  viscous 
dissipation  and  the  potential  energy  of  the  upper  ocean  remains  constant. 

Heat  dominated  mixing  is  characterized  by  a positive  heat  flux  at 
the  ocean's  surface.  The  depth  of  the  mixed  layer  is  determined  primar- 
ily by  the  magnitude  of  the  surface  heat  flux  and  the  strength  of  the 
wind  stress.  If  there  is  a positive  heat  flux  at  the  sea  surface,  the 
wind  stress  can  do  work  on  the  upper  ocean  by  mixing  this  heat  through- 
out the  mixed  layer  as  well  as  by  deepening  the  mixed  layer  and  mixing 
dense  water  upwards.  Since  continued  deepening  of  the  mixed  layer  is  not 
necessary  in  order  for  the  wind  stress  to  do  work  on  the  upper  ocean, 
the  Denman  and  Kraus  and  Turner  models  as  well  as  the  three  models  com- 
pared here  predict  that  a stationary  mixed  layer  depth  is  established  if 
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the  wind  stress  and  surface  heat  flux  remain  fairly  constant  in  time. 

Because  of  the  simple  nature  of  the  equations  and  the  boundary  and 
initial  conditions  for  both  problems  of  mixed  layer  deepening,  the  prob- 
lems are  amenable  to  a similarity  scaling  which  results  in  a set  of  equa- 
tions involving  only  one  free  parameter.  Hence,  in  both  problems  the 
fairly  stationary  mixed  layer  depths  achieved  after  several  days  of  inte- 
gration are  a function  of  a single  nondimensional  parameter  for  each  of 
the  models.  This  considerably  simplifies  the  analysis  of  the  problems 
and  provides  some  insight,  not  only  into  how  the  mixed  layer  depth  de- 
pends on  the  surface  wind  stress  and  heat  flux,  the  stratification,  and 
the  local  inertial  frequency,  but  also  into  why  the  mixed  layer  depths 
predicted  by  the  diffusion  models  differ.  Mellor  and  Durbin  (197*0 
performed  such  a scaling  for  the  problem  of  the  stationary  mixed  layer 
depth  determined  by  a constant  positive  surface  heat  flux,  and  noted  that 
the  mixed  layer  depth  appeared  to  be  given  by 


2 !3gQ 


(1) 


in  the  limit  U - Here,  p is  the  ratio  of  the  turbulent  Ekman  layer 
length  scale  u /f  and  the  Monin-Obukoff  length  scale  uVBgQ,  u is  the 
friction  velocity  associated  with  the  wind  stress,  0 is  the  coefficient 
of  thermal  expansion  of  sea  water,  g is  the  gravitational  constant,  and 
Q is  the  surface  heat  flux.  It  should  be  noted  that  the  mixed  layer 
depth  does  not  depend  upon  the  Coriolis  parameter  f in  this  limiting  case 
Models  of  the  mixed  layer  fall  into  two  general  groups,  integral 
models  and  grid  point  models.  Integral  models  such  as  those  of  Kraus  and 
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Turner  (1967)  and  Denman  (1973)  treat  the  mixed  layer  as  a vertically 
integrated  layer  and  use  the  vertically  averaged  values  of  mixed  layer 
quantities  to  predict  the  behavior  of  the  mixed  layer.  Grid  point  models 
use  finite  differenced  forms  of  the  conservation  equations  over  a spa- 
tially fixed  grid  to  predict  the  vertical  distributions  of  velocity  and 
density  with  time.  These  models  have  much  in  common  with  the  grid  point 
atmospheric  boundary  layer  models,  and  some  models,  such  as  those  of  Vager 
and  Zilitinkevich  (1968)  and  Mellor  and  Yamada  (1975)  have  been  used  to 
model  both  atmospheric  and  oceanic  boundary  layers. 

Integral  models  use  depth  averaged  values  for  the  density  and 
velocity  of  the  mixed  layer,  and  this  is  justified  on  the  grounds  that 
turbulent  mixing  within  the  mixed  layer  rapidly  mixes  the  layer  to  a 
fairly  uniform  composition.  However,  although  the  rate  of  mixing  within 
the  mixed  layer  is  rapid,  it  is  not  instantaneous.  For  instance,  if  the 
depth  of  the  mixed  layer  is  about  30  meters  and  the  effective  turbulent 
diffusivity  is  about  100  cm2/sec  (Ostapoff  and  Worthem,  1972),  the 
diffusive  time  scale  within  the  mixed  layer  is  about  a day,  and  it  will 
be  seen  that  the  diffusion  models  compared  here  predict  that  the  mixed 
layer  depth  cap  take  a couple  of  days  or  more  to  adjust  to  a change  in 
the  sea  surface  heat  flux.  The  temperature  and  velocity  structure  with- 
in the  mixed  layer  affects  the  response  of  the  mixed  layer  to  changes  in 
the  wind  stress  and  heat  flux  at  the  surface. 

Grid  point  models  predict  the  vertical  distribution  of  velocity 
and  density  almost  continuously  in  space.  Hence,  they  model  the  diffu- 
sion of  heat  and  momentum  within  the  mixed  layer  as  well  as  the  fluxes  of 
heat  and  momentum  at  the  surface  and  at  the  bottom  of  the  mixed  layer  and 
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below.  Another  advantage  of  grid  point  models  Is  that  no  assumptions 
need  to  be  made  concerning  the  behavior  of  the  mixed  layer  Itself.  The 
turbulent  diffusivity  is  determined  by  the  mean  velocity  and  density 
fields  and  the  vertical  distribution  of  the  turbulence  in  turn  defines 
the  extent  of  the  mixed  layer. 

The  advantage  of  Integral  mixed  layer  models  is  their  simplicity 
and  speed  of  computation.  The  limitations  of  computation  time  in  three- 
dimensional  numerical  simulations  of  the  ocean  may  prohibit  the  use  of 
iSixed  layer  models  requiring  high  vertical  resolution.  However,  for  one- 
dimensional  simulations  of  the  upper  ocean  where  horizontal  inhomogene- 
ities are  Ignored,  the  computation  time  required  for  integrations  of  days 
or  months  with  any  but  the  most  complicated  grid  point  models  is  not  sub- 
stantial, and  these  models  are  to  be  preferred  over  integral  models. 

Increased  understanding  of  the  dynamics  of  the  mixed  layer  through 
observation  and  simulation  with  grid  point  models  may  allow  the  opti- 
mization of  Integral  mixed  layer  models  for  particular  oceanic  applica- 
tions. 

II.  DESCRIPTION  OF  THE  UPPER  OCEAN 

The  temperature  structure  of  the  upper  ocean  has  a major  effect  on 
the  vertical  turbulent  diffusivity.  The  upper  ocean  is  warmed  by  solar 
radiation,  exchanges  heat  with  the  atmosphere  through  sensible  heat  ex- 
change, and  loses  heat  due  to  short  wave  back  radiation  and  evaporation. 
During  the  spring  and  summer  heating  season,  solar  heating  predominates 
and  a stable  density  gradient  is  built  up  in  the  upper  layers.  Winds 
blowing  over  the  sea  generate  turbulent  mixing  in  the  near  surface  water. 


forming  a layer  of  fairly  uniform  temperature  and  salinity.  The  mixing 
of  the  near  surface  water  causes  a large,  stable  temperature  gradient, 
a thermocline,  to  develop  between  the  wind  mixed  layer  and  the  water  be- 
low. This  temperature  gradient  inhibits  vertical  mixing  and  tends  to 
prevent  turbulence  from  penetrating  below  the  thermocline.  Deepening 
of  the  mixed  layer  due  to  increasing  winds  or  wanting  of  the  layer  through 
heating  at  the  surface  act  to  strengthen  the  thermocline,  creating  an  even 
stronger  barrier  to  the  penetration  of  heat  and  momeutun  into  deeper  water. 

The  effective  vertical  diffusivity  in  the  upper  ocean  can  vary  from 
values  on  t*-  order  of  100  cm^/sec  in  the  turbulent,  wind  mixed  surface 
layer  to  values  approaching  the  molecular  rates  of  diffusion  of  heat  and 
momentum  in  nonturbulant  regions  below  the  thermocline. 

III.  DESCRIPTION  OF  THREE  EDDY  DIFFUSION  MODELS 

The  first  significant  modeling  of  a diffusive  boundary  layer  in 
the  presence  of  rotation  was  by  Ekman  (1905).  Ekman  analysed  the  re- 
sponse of  a homogeneous  ocean  with  constant  vertical  diffusion  coefficient 
K to  a suddenly  imposed  surface  wind  stress  t.  He  found  that  because  of 
the  earth's  rotation  the  frictional  influence  of  the  wind  was  limited  to 
a boundary  layer  of  depth  rr'|  2 K/f,  where  f is  the  Coriolis  parameter. 

This,  the  Ekman  layer  depth,  is  the  depth  at  which  the  current  has  been 
reduced  to  e ^ times  its  value  at  the  surface.  Ekman  noted,  however, 
that  the  rate  of  vertical  diffusion  in  the  real  ocean  would  depend  upon 
the  stratification. 
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Munk- Anderson  Diffusion  Modal 


Munk  end  Anderson  ( 19^3)  reasoned  from  dimensional  arguments 


that  the  effective  vertical  diffuaivity  In  a turbulent,  stratified  fluid 


should  depend  strongly  upon  the  Richardson  number  Ri,  where 


Here  g is  the  acceleration  of  gravity  o is  the  potential  density  of  the 
fluid,  z is  the  vertical  coordinate,  and  u and  v are  the  horizontal 


density  stratification. 

Munk  and  Anderson  postulated  the  functional  forms 


for  the  vertical  diffusion  coefficients  for  momentum  and  heat  and 

Here  K is  the  magnitude  of  the  vertical  diffusion  coefficients  in  the 
o 

limit  Ri  - 0 and  r^,  i^,  &M,  and  0H  are  constants.  From  an  analysis  of 
wind  profile  measurements  in  the  atmospheric  boundary  layer,  Munk  and 
Anderson  determined  that  these  constants  should  be  assigned  the  values 
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B.  Vager-Zilitinkevlch  Diffusion  Model 

Vager  and  Zilitinkevich  ( 19^8)  developed  a diffusion  model 
based  upon  Che  assumption  that  the  rate  of  vertical  diffusion  of  heat  and 
momentum  is  proportional  to  the  square  root  of  the  kinetic  energy  of  the 
turbulent  eddies  B and  to  the  length  scale  of  the  turbulence  L,  so  that 

>>,  - "h  - - W 


where  cq  is  a constant.  In  their  model,  the  kinetic  energy  of  the  turbu- 
lence is  predicted  by  a simplified  form  of  the  equation  for  the  turbulent 
kinetic  energy 

M.oA/v  + :Vgs.c£, 

at  ~ at  \ " at  / * [\  at  / \ at  / J “ at  k. 


where  c is  a universal  constant  with  value  0.04-6,  0£g  is  a constant  with 
value  0.73,  and  0 is  the  coefficient  of  thermal  expansion  of  the  fluid. 
The  terms  on  the  right  side  of  Eq.  (5)  model, respectively,  the  diffusion 
and  production  of  turbulent  energy,  the  conversion  of  turbulent  energy  to 
potential  energy,  and  the  dissipation  of  turbulent  energy.  The  constants 
c and  cq  are  related  by  cq  ■ c1^4.  In  solving  Eq.  (5),  the  flux  of 
turbulent  kinetic  energy  at  the  ocean's  surface  is  taken  to  be  zero,  so 
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0 at  z ■ 0 . 
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The  length  scale  l is  determined  from  the  hypothesis  of  Zilitinkevich  and 
Laykhtman  concerning  the  length  scale  of  turbulence  In  a turbulent  bounda- 
ry layer,  in  which 


l • - k 


±H' 

dz  l 


where  k is  von  Groan's  constant,  taken  to  be  .h.  This  expression  can  be 
integrated  to  the  more  convenient  form 


I-*-  k + f f . 
K J >Im 


where  z is  the  roughness  parameter  typical  of  the  boundary  and  B is  the 
o ° 

turbulent  kinetic  energy  at  the  boundary.  Combining  Eqs.  (h)  and  (8), 
the  vertical  diffusion  coefficients  can  be  expressed  as 
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C.  Mellor-Yamada  Level  2 Diffusion  Model 

Mellor  and  Yamada  (197*0  have  developed  a series  of  four 
turbulence  closure  models  for  planetary  boundary  layers  labeled  Level  1 
to  Level  4 in  order  of  increasing  complexity.  The  models  were  developed 
through  systematic  scaling  of  the  equations  for  the  turbulent  fluxes  of 
heat  and  momentum.  The  modeling  of  the  triple  correlations  and  dissipa- 
tion terms  in  these  equations  was  based  on  hypotheses  proposed  by  Rotta 
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and  Kolmogoroff  and  Che  values  of  cons cants  used  in  such  terms  were  de- 
termined from  unstratified  turbulence  data. 

Because  of  its  simple  form,  the  Level  2 model  was  chosen  for  com- 
parison with  the  Munk-Anderson  and  Vager-Zilitinkevich  models.  For  the 
Level  2 model,  the  turbulence  fluxes  of  momentum  and  heat  can  be  ex- 
pressed in  eddy  coefficient  form  where  the  eddy  coefficients  depend  upon 
the  local  mean  velocity  and  density  gradients.  Henceforth,  reference  to 
the  Mellor-Yamada  model  will  refer  to  their  Level  2 turbulence  closure 
model.  Melloz  and  Durbin  (1975)  have  applied  this  model  to  the  mixed 
layer  of  the  ocean. 

The  vertical  diffusion  coefficients  for  the  Level  2 model  are  de- 
fined as 


and 


where  the  stability  functions  and  are  given  by 


SH  - 0.55664  - 1.97808  Rf/(l  - Rf)  , 


(10) 


and 


(ID 


Sh 


(0.52  - 1.404  R /(I  - R )) 

s - -I  — -.  - z ■ - - 

H (0.688  - 2.068  R./(l  - R,)) 


The  flux  Richardson  number  Rf  is  related  to  the  Richardson  number  (2)  by 
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The  turbulent  energy  B is  calculated  from  a highly  simplified  turbulent 
energy  equation  given  by 
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Equation  (13)  expresses  a balance  among  the  production  of  turbulent  energy 
from  the  mean  velocity  shear,  the  conversion  of  turbulent  energy  to  po- 
tential energy,  and  the  dissipation  of  turbulent  energy.  The  turbulence 
length  scale  L is  determined  from  the  vertical  extent  of  the  turbulence 
field  as 
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The  differences  in  the  three  diffusion  models  can  be  seen  most 
readily  by  considering  how  the  turbulent  diffusion  coefficients  depend 
upon  the  Richardson  number. 

For  the  Munk- Anders on  model , 
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Hence , 
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where  and  IC,  are  the  diffusion  coefficients  for  an  unstratified 
o o 

ocean  where  Ri  is  zero. 


For  the  Vager-Zilitinkevich  model 
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Hence  , 
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For  the  purpose  of  deriving  Eq.  (17)  the  terms  representing  the  storage 
and  diffusion  of  turbulent  kinetic  energy  in  Eq.  (5)  were  ignored.  In 
actual  integrations  using  the  Vager-Zilitinkevich  model,  these  terms  were 
not  found  to  affect  predicted  mixed  layer  depths  significantly. 

For  the  Mel lor  Level  2 model 


and 
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Hence 
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Equations  (16),  (19)»  and  (21)  are  plotted  in  Fig.  1.  The  critical 
Richardson  number  Ri.  for  which  turbulent  mixing  ceases  differs  signifi- 
cantly for  each  of  the  three  models.  The  Vager-Zilitinkevich  and  Mellor- 
Yamada  models  predict  that  turbulence  will  be  "shut'off"  for  Richardson 
numbers  exceeding  1.0  and  .23  respectively.  In  contrast,  the  Munk- 
Anderson  model  predicts  a very  gradual  decrease  in  the  intensity  of  the 
turbulent  mixing  as  the  Richardson  number  increases. 

I V.  ADJUSTMENT  OF  THE  DIFFUSION  MODELS 

For  a given  set  of  boundary  and  initial  conditions  it  was  felt  that 
the  three  diffusion  models  should  predict  vertical  diffusivities  within 
the  mixed  layer  of  about  the  same  magnitude  as  the  vertical  diffusivi- 
ties that  have  been  observed  in  the  upper  ocean  under  similar  conditions. 
Bowden,  Howe,  and  Tait  (1970)  derived  estimates  of  the  diffusivity  of 
heat  Kjj  of  about  93  cm^/sec  in  the  near  surface  layer  from  diurnal 
temperature  variations,  and  Ostapoff  and  Northern  (1972)  deduced  an 
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(subscript  M)  and  heat  (subscript  H)  for  the  three  diffusion  models 


average  value  for  of  about  72  cs^/sec  in  the  mixed  layer  during  a day 
when  the  water  was  being  stably  stratified  by  solar  heating  and  the  wind 
speed  was  in  the  range  of  6 to  7 m/sec.  Since  the  conditions  under  which 
these  observations  were  made  were  similar  to  the  conditions  imposed  on 
the  diffusion  models  for  the  two  problems  of  mixed  layer  deepening,  where 
the  wind'stress  was  one  dyne  and  the  mixed  layer  depth  was  fairly  shallow, 
about  30  meters,  the  diffusion  models  were  adjusted  to  predict  a maximum 
diffusivity  of  about  100  cn^/sec. 

Munk  and  Anderson  ( 19^8)  proposed  that  the  parameter  Kq  in  their 
model  (Eq.  2)  should  be  a function  only  of  the  wind  stress,  and  their 
paper  gives  a value  for  Kq  of  150  caP/sec  for  a wind  stress  of  one  dyne. 
In  the  upper  ocean,  it  might  be  expected  that  Kq  should  depend  upon  the 
depth  of  the  mixed  layer  as  well  as  the  wind  stress  since  as  the  mixed 
layer  deepens,  the  length  scale  of  the  turbulent  eddies  and  the  magnitude 
of  the  vertical  diffusivity  would  be  expected  to  Increase.  Hence,  for  a 
given  wind  stress,  a particular  value  for  Kq  might  be  too  small  for  a 
very  deep  mixed  layer  and  too  large  for  a very  shallow  mixed  layer.  For 
the  two  problems  of  mixed  layer  deepening  studied  here,  Kq  was  set  equal 
to  100  cm2 /sec. 

The  Vager-Zilitinkevich  model  as  described  in  Section  3 predicted 
vertical  diffusiviti es  greater  than  300  anP/sec  for  the  problems  studied 
which  seemed  too  large,  so  the  turbulence  length  scale  as  predicted  by 
Eq.  (7)  was  reduced  by  a factor  of  three.  The  Mellor-Yamada  model  pre- 
dicted diffusivities  in  fair  agreement  with  the  observed  values  and  no 
adjustment  to  it  was  made. 
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It  must  be  noted  that  the  adjustments  to  the  magnitudes  of  the 
eddy  diffusivities,  which  did  not  affect  the  functional  dependence  of 
the  diffusivities  upon  the  stratification,  did  not  significantly  change 
the  mixed  layer  depths  predicted  by  the  diffusion  models  for  the  two 
mixed  layer  deepening  problems  studied. 

Hell  below  the  mixed  layer  the  diffusion  models  predict  no  turbu- 
lent mixing,  and  thus  the  molecular  diffusivities  assert  themselves. 
Although  observations  indicate  that  intermittent  turbulent  mixing  en- 
hances vertical  diffusion  below  the  mixed  layer  (Grant,  Moillet,  and 
Vogel,  1968),  mechanisms  for  the  generation  of  turbulence  in  this  region 
(such  as  the  breaking  of  internal  waves)  were  ignored  here.  The  mixed 
layer  depths  predicted  by  the  models  were  not  found  to  be  sensitive  to 
the  eddy  diffusivity  imposed  below  the  mixed  layer  provided  it  was  small 
relative  to  the  rate  of  mixing  within  the  mixed  layer. 

V.  ONE-DIMENSIONAL  UPPER  OCEAN  MODEL  EQUATIONS 

The  equations  for  the  conservation  of  heat  and  momentum  for  a one- 
dimensional model  of  the  upper  ocean  where  the  only  spatial  variations 
are  with  respect  to  the  vertical  coordinate  z can  be  written  as 


A small  up welling  velocity  typical  of  the  open  ocean  could  have  been  in- 
cluded but  was  not  considered  to  be  important  to  the  dynamics  of  the 
mixed  layer  on  the  time  scale  of  the  problems  considered  here.  The 
salinity  of  the  upper  ocean  was  assumed  to  be  uniform  and  the  density 
of  the  seawater  was  taken  to  be  linearly  dependent  upon  the  temperature, 
the  coefficient  of  thermal  expansion  of  seawater  & being  2.3  x 10”4  °C_1 . 

The  boundary  conditions  for  Eqs.  (22)  can  be  written  as 
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where  u^  ■ (r/p)1^2  is  the  friction  velocity  associated  with  the  wind 
stress  t,  and  Q is  the  sea  surface  heat  flux. 

For  the  first  problem,  the  case  of  a constant  imposed  wind  stress, 
zero  surface  heat  flux,  and  a linear  initial  stratification,  u^  was  taken 
to  be  1 cm/sec  and  Q was  taken  to  be  zero,  and  at  t ■ 0 the  ocean  was 
assumed  to  be  at  rest  and  stratified  with  a linear  temperature  gradient 
of  .1  °C/meter. 

For  the  second  problem,  the  case  of  a constant  imposed  wind 
stress  and  a constant  positive  surface  heat  flux,  u^  was  taken  to  be 
1 cm/sec  and  Q was  taken  to  be  3*55  x 10  C-cm/sec  which  corresponds 
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to  a surface  heat  flux  of  300  cal/cm2/day  into  the  ocean.  At  t - 0 the 
ocean  was  taken  to  be  at  rest  and  unstratified. 

VI.  SCHEME  FOR  NUMERICAL  SOLUTION  OF  THE  OCEAN  MODEL  EQUATIONS 

The  Eqs.  (22)  for  the  mean  velocity  and  temperature  fields  were 
solved  numerically  by  finite  differencing  in  space  and  time.  A 2- level 
time  differencing  scheme  was  used.  The  diffusion  terms  were  evaluated 
at  the  n + 1 time  level  for  stability  and  the  Coriolis  terms  were  split 
equally  between  the  n and  n + 1 time  levels.  The  diffusion  coefficients 
were  evaluated  at  time  level  n. 

A uniform  spatial  grid  was  used  with  a spacing  of  two  meters 
between  grid  points.  The  temperature,  density,  and  velocity  were  speci- 
fied at  each  grid  point  and  the  eddy  coefficients  were  calculated  at 
positions  between  the  grid  points.  A time  step  of  six  minutes  was  used 
for  the  integrations. 

VII.  RESULTS  OF  NUMERICAL  INTEGRATION  OF  THE  DIFFUSION  MODELS 

A.  Wind  Dominated  Deepening  of  the  Mixed  Layer 

Figures  2 to  6 summarize  the  results  of  the  integration  of 
the  three  diffusion  models  for  the  first  problem,  the  case  of  a constant 
wind  stress  of  1 dyne/cm2,  zero  surface  heat  flux,  and  an  initial  strati- 
fication of  .1  °C/meter. 

The  deepening  of  the  mixed  layer  with  time  is  shown  in  Fig.  2. 

For  the  Vager-Zilitinkevich  and  Mellor-Yamada  models,  the  mixed  layer  is 
fairly  well  defined,  since  the  production  of  turbulent  energy  by  the  mean 
shear  stops  abruptly  when  the  Richardson  number  exceeds  its  critical 
value  Ric*  the  Munk-Anderson  model,  the  turbulent  mixing  is  never 
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Fig.  2 — The  mixsd  layer  depth  predicted  by  the  three  diffusion  models  for 
the  problem  of  wind  dominated  deepening.  The  initial  stratification  was  0.1 
°C/meter.  The  imposed  wind  stress  was  1 dyne/cm2  and  the  surface  heat 
flux  was  set  to  zero. 
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completely  shut  off  by  the  stratification.  Hence,  for  the  purpose  of  cot^ 
structing  Fig.  2,  the  mixed  layer  was  taken  to  be  the  region  of  the  upper 
ocean  for  which  the  turbulent  diffusivity  for  heat  was  greater  than 
1.0  cs^/sec.  Figure  2 shows  that  the  downward  mixing  proceeded  rapidly 
at  first,  but  slowed  considerably  after  about  one  inertial  period  (36 
hours).  Although  a truly  steady  state  was  never  reached,  condi- 
tions within  the  mixed  layer  became  nearly  stationary  after  a few  days 
and  further  erosion  of  the  thermocline  proceeded  very  slowly. 

Figure  3 displays  the  variation  of  the  vertical  diffusion  coef- 
ficients with  depth  after  the  models  had  run  for  13  days.  The  Vager- 
Zilitinkevich  and  Mellor-Yamada  models  effectively  shut  off  turbulent 
mixing  at  the  thermocline.  The  Munk-Anderson  model  does  not,  and  there- 

to 

fore  their  model  predicted  continued  deepening  of  the  mixed  layer  even 
after  several  days. 

The  temperature  profiles  in  Fig.  ^ show  the  formation  of  the 
thermocline  and  the  mixed  layer.  After  several  days  of  mixing  with  no 
surface  heat  flux  the  mixed  layer  was  fairly  uniform  in  temperature. 

The  Vager-Zilitinkevich  and  Mellor-Yamada  models  show  a much  sharper 
thermocline  than  the  Munk-Anderson  model  where  the  thermocline  was  still 
being  eroded  by  turbulence  after  13  days. 

The  amplitude  of  the  inertial  oscillations  within  the  mixed  layer 
generated  by  the  sudden  change  in  the  wind  stress  at  t ■ 0 is  shown  in 
Fig.  3*  The  inertial  oscillations  within  the  mixed  layer  are  constant 
in  amplitude  and  phase,  so  that  the  entire  mixed  layer  oscillates  as  a 
'slab'  over  the  thermocline.  These  Inertial  oscillations  persisted  with 
little  change  in  amplitude  over  the  13  days  of  integration  because  of 
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the  very  small  (and  probably  unrealistic)  damping  due  to  the  small  value 
of  below  the  mixed  layer.  The  amplitude  of  Inertial  oscillations  in 
the  upper  ocean,  besides  being  dependent  upon  changes  in  the  direction 
and  strength  of  the  wind  (Pollard  and  Millard,  1970),  is  roughly  in- 
versely proportional  to  the  mixed  layer  depth  and  this  relationship  can 
be  seen  by  comparing  the  mixed  layer  depths  and  inertial  oscillation 
amplitudes  predicted  by  the  three  models. 

The  mean  velocity  profiles  are  shown  in  Fig.  6.  These  were  the 
mean  velocities  at  each  depth  about  which  the  inertial  oscillations 
occurred.  A velocity  profile  for  a classical,  steady,  homogeneous  Ekman 
layer  with  » 100  cn^/sec  is  included  in  Fig.  6 for  comparison.  An 
analysis  of  Eqs.(22)  shows  that  the  Ekman  layer  transport  depends  only 

on  the  wind  stress  and  the  Coriolis  parameter  and  is  directed  90°  to  the 

* 

right  of  the  wind  stress.  Therefore,  all  the  models  yie^l  the  same 
transport  as  that  predicted  by  Ekman' s (1905)  theory,  in  which 
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B.  Heat  Dominated  Deepening  of  the  Mixed  Layer 


Figure  7 shows  the  change  of  the  mixed  layer  depth  with  time 
predicted  by  the  three  diffusion  models  for  the  second  problem,  the  case 
of  an  imposed  wind  stress  of  1 dyne/cm2,  an  imposed  surface  heat  flux 
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of  300  cal/cm2/day,  and  negligible  initial  stratification.  The  qualita- 
tive behavior  predicted  by  the  three  models  was  similar.  The  mixed  layer 
mixed  very  deeply  at  first  because  of  the  absence  of  stratification,  but 
as  the  stratification  of  the  upper  ocean  increased  due  to  the  downward 
mixing  of  the  heat  entering  at  the  surface,  the  deepening  slowed  and 
stopped.  As  the  stratification  continued  to  increase,  the  turbulent 
mixing  was  weakened  still  further,  and  the  mixed  layer  depth  decreased 
until  a certain  balance  was  reached  between  the  production  of  turbulent 
energy  by  the  mean  shear  and  the  production  of  potential  energy  by  the 
turbulent  mixing.  At  this  point  the  mixed  layer  depth  stabilized  since 
the  density  gradient  was  just  sufficient  to  stabilize  the  mean  velocity 
shear  at  the  base  of  the  mixed  layer. 

Because  of  the  larger  value  of  the  critical  Richardson  number  of 
the  Vager-Zilitinkevich  model,  this  model  predicted  a deeper  mixed  layer 
than  the  Mellor-Yamada  model  for  both  problems  of  mixed  layer  deepening. 

A larger  value  of  the  critical  Richardson  number  implies  that  a stronger 
density  gradient  is  required  to  stabilize  a given  velocity  shear.  As  the 
mixed  layer  deepens,  the  mean  velocity  shear  near  the  base  of  the  mixed 
layer  decreases.  Hence,  the  larger  the  value  of  the  critical  Richardson 
number,  the  deeper  the  mixed  layer  can  penetrate  before  the  turbulent 
mixing  at  the  base  of  the  mixed  layer  becomes  suppressed  by  the  strati- 
fication. 

Figure  7 shows  that  the  Munk-Anderson  model  took  a longer  time  to 
establish  a steady  mixed  layer  depth  for  the  second  problem  of  mixed 
layer  deepening  than  the  two  other  diffusion  models.  The  mixed  layer 
depth  for  the  Munk-Anderson  model  was  not  very  well  defined  until  about 
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Fig.  7 — The  mixed  layer  depth  predicted  by  the  three  diffusion  models  for  the 
problem  of  heat  dominated  mixing.  The  upper  ocean  was  initially  unstratified. 
The  imposed  wind  stress  was  1 dyne/cm2  and  the  imposed  sea  surface  heat  flux 
was  300  cal/cm2/day. 


Che  Centh  day  of  integration,  and  so  the  mixed  layer  depth  up  to  that 
time  is  shown  with  a dotted  line.  A steady  mixed  layer  depth  becomes 
established  with  the  Munk-Anderson  model  for  problems  of  heat  dominated 
deepening  because,  even  though  some  mixing  occurs  below  the  mixed  layer, 
there  is  a certain  depth  at  which  the  downward  transport  of  heat  cannot 
keep  pace  with  the  heat  flux  at  the  surface.  Hence,  the  stratification 
builds  up  at  this  depth,  and  mixing  below  is  gradually  shut  off  as  the 
turbulence  finds  it  increasingly  difficult  to  penetrate  the  growing 
stratification. 

VIII.SIMILAIi.ITY  SCALING  ANALYSIS  OF  TWO  PROBLEMS  OF  MIXED  LAYER  DEEPENING 

A similarity  scaling  analysis  of  the  problems  of  wind  dominated 
and  heat  dominated  deepening  previously  described  can  provide  some  in- 
sight into  the  dependence  of  the  predicted  mixed  layer  depth  upon  the 
problem  parameters.  Because  they  can  be  nondimens ionalized  in  the  same 
fashion,  the  Vager-Zilitinkevich  and  Mellor-Yamada  models  will  be  in- 
cluded in  this  analysis. 

Since  the  two  problems  of  mixed  layer  deepening  were  defined  very 
simply,  they  can  be  scaled  so  that  the  only  variable  parameter  remaining 
in  either  of  the  problems  occurs  in  the  formulation  of  the  eddy  coeffi- 
cients of  heat  and  momentum  due  to  the  scaling  of  the  Richardson  number. 
For  the  problem  of  wind  dominated  deepening  with  a constant  wind  stress, 
no  surface  heating,  and  a linear  initial  stratification,  this  parameter 
is 
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which  is  Che  square  of  Che  raCio  of  Che  Brunt-Vaisala  frequency  N and 
Che  local  inerCial  frequency  f.  Here,  N "\jfSgr,  where  T is  che  initial 
stratification  of  Che  upper  ocean,  0 is  Che  coefficient  of  thermal  ex- 
pansion of  che  seawater;  and  g is  the  gravitational  constant.  For  the 
problem  of  heat  dominated  deepening  with  a constant  wind  stress,  a con- 
stant positive  heat  flux  Q,  and -negligible  initial  stratification,  the 
variable  parameter  remaining  after  scaling  is 

0'  (26) 


which  is  the  raeio  of  che  turbulent  Ekman  layer  length  scale  u / f and 
the  Monin-Obukov  length  scale  u^_/0gQ,  where  u^  is  the  friction  velocity 
associated  with  che  wind  stress. 

This  scaling  can  be  made  more  general  if  approximations  are  made 
to  the  Vager-Zilitinkevich  and  Mellor-Yamada  diffusion  models  to  allow 
two  additional  parameters  Co  be  brought  into  che  scaling  which  are  charac- 
teristic of  the  diffusion  models  themselves.  One  of  these  parameters  is 
a,  a constant  to  which  the  eddy  coefficients  are  proportional  and  which 
can  be  adjusted  to  change  the  magnitude  of  the  eddy  coefficients  without 
affecting  the  nature  of  the  dependence  of  the  eddy  coefficients  upon  the 
stratification.  The  other  parameter  is  Ri£,  which  is  the  critical  Richard- 
son number  for  each  of  the  models  above  which  the  generation  of  turbu- 
lence by  the  mean  shear  is  completely  suppressed  by  the  stratification. 

Although  a and  Ri£  are  fixed  constants  for  both  the  Vager- 
Zilitinkevich  and  Mellor-Yamada  diffusion  models,  by  treating  them  as 


having  a degree  of  variability  some  assessment  can  be  made  of  the 
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Importance  of  these  constants  in  determining  the  mixed  layer  depths  pre- 
dicted by  the  two  models. 

The  approximations  are  made  to  the  models  in  order  to  simplify  the 
scaling  and  at  the  same  time  present  the  models  in  a form  that  makes 
their  similarities  and  differences  more  apparent.  Both  the  Vager- 
Zilitinkevich  and  Mellor-Yamada  models  can  be  fairly  well  approximated 
by  the  general  expression 

Kjj  - Kjj  - 0UZ  |vj  (1  - Ri/Ric)Y  , (27) 

where  and  are  the  eddy  diffusivities  for  heat  and  momentum,'  & is 

the  turbulence  length  scale,  and  |Vz|  is  the  absolute  value  of  the  mean 

velocity  shear.  For  the  Vager-Zilitinkevich  model  y ■ 1/2,  a ■ 1.0,  and 

Ri  * 1.0,  and  for  the  Mellor-Yamada  model  y =•  3/2,  a 1.0,  and  Ri  * 
c c 

0.23.. 

Putting  the  Vager-Zilitinkevich  and  Mellor-Yamada  models  into  the 
form  of  Eq.  (27)  requires  dropping  the  tendency  and  diffusion  terms  in 
the  turbulent  energy  Eq.  (5)  for  the  Vager-Zilitinkevich  model  and 
simplifying  the  functions  S^(Ri)  and  S^(Ri)  of  the  Mellor-Yamada  model 
Eq.  (ll)  to  the  form  (l  - Ri/Ric)Y.  The  results  of  a test  integration 
of  the  approximate  form  of  the  models  differed  little  from  the  results 
of  integration  of  the  actual  models.  It  will  be  mentioned  here  that  all 
the  results  presented  were  obtained  from  integration  of  the  diffusion 
models  in  their  proper  form  as  described  in  Sections  III  and  IV. 

For  an  unstratified  flow  where  the  Richardson  number  is  zero, 
both  the  Vager-Zilitinkevich  and  Mellor-Yamada  models  essentually  reduce 
to  Prandtl's  mixing  length  theory  in  which  the  eddy  coefficients  are 
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proportional  to  the  square  of  the  turbulence  length  scale  and  the  mean 
velocity  shear.  For  a stratified  flow  the  eddy  coefficients  are  modified 
by  a function  of  the  Richardson  number  which  reduces  the  vertical  dif- 
fusivity  for  a given  velocity  shear  as  the  Richardson  number  increases, 
and  hence,  models  the  suppression  of  the  turbulence  by  the  stratification. 
For  the  approximate  form  of  the  Vager-Zilitinkevich  and  Mellor-Yamada 
models  given  by  Eq.  (27),  this  stratification  function  has  the  form 


(28) 


The  values  of  the  constants  Ri£  and  y in  Eq.  (28)  differentiate 
the  Vager-Zilitinkevich  and  Mellor-Yamada  models.  The  most  significant 
difference  between  the  models  is  that  the  critical  Richardson  number  for 
the  Vager-Zilitinkevich  model  Ric  * 1.0  is  much  higher  than  that  for  the 

Mellor-Yamada  model  for  which  Ric  ■ .23.  A probably  less  significant 
difference  is  that  the  ’shape'  of  the  stratification  function  (28)  which 
is  determined  by  the  value  of  y is  quite  different  for  the  two  models. 

The  Vager-Zilitinkevich  model  with  y « 1/2  predicts  considerably  less 
suppression  of  the  turbulent  mixing  for  a given  value  of  Ri/Ri^  than 
does  the  Mellor-Yamada  model  with  y *•  3/2.  Since  the  parameter  Ric  is 
to  be  included  in  the  scaling  of  the  problem  of  mixed  layer  deepening, 
differences  in  the  nondimensionalized  mixed  layer  depth  predicted  by  the 
models  should  be  due  to  the  difference  in  the  'shape'  of  the  stratifica- 
tion function  (28)  represented  by  the  different  values  of  y. 

The  scaling  used  to  nondimens ionalize  the  two  problems  of  mixed 
layer  deepening  is  similar  to  that  suggested  by  Mellor  and  Durbin  (197*0. 
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A . Wind  Dominated  Deepening 

For  the  case  of  wind  dominated  deepening  with  a constant 
wind  stress,  zero  surface  heat  flux,  and  a linear  initial  temperature 
stratification  T,  the  depth  can  be  scaled  as  u^^cT/f,  the  time  as  1/f, 
the  velocity  as  u^/  \|^  the  temperature  as  ru^Nfa/f,  and  the  eddy  coeffi- 
cients as  u^y/a/f.  After  scaling,  the  only  parameter  remaining  in  the 
problem  is 


- jfsl 

1 f2Ri 


0 


(29) 


where  N “ ^&gr  is  the  Brunt-Vaisala  frequency  determined  by  the  initial 
stratification  of  the  upper  ocean.  Therefore,  the  fairly  steady  mixed 
layer  depth  h achieved  after  several  days  of  integration  is  a function 


of  0^  and  is  given  by 


uja 


h*(  0 ) , 


(30) 


where  h#  is  the  dimensionless  mixed  layer  depth. 

Figure  8 shows  the  dependence  of  h*  upon  0 found  from  numerical 
integration  of  the  Mellor-Yamada  and  Vager-Zilitinkevich  models.  The 
two  curves  almost  coincide, indicating  that  most  of  the  difference  in 
mixed  layer  depth  predicted  by  the  two  models  for  the  case  of  wind 
dominated  deepening  is  due  to  the  difference  in  Ri£  (which  is  taken  into 
account  by  the  scaling).  The  difference  in  the  value  of  y in  the  strati- 
fication function  (26)  of  the  two  models  and  the  diffusion  of  turbulent 
energy  in  the  Vager-Zilitinkevich  model  did  not  significantly  affect  the 
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predicted  mixed  layer  depths. 


For  values  of  0^  greater  than  105  (which  corresponds  to  stratifi- 
cations exceeding  10  3 °C/an  when  the  Coriolis  parameter  is  about  5 x 10"s 
sec  1),  the  slope  of  h#(0^)  is  about  - 1/4  and  h*  is  approximately  given 
by 


-i/4 

h*  =-  2.7  0 


(5D 


Therefore, 


h 


2.7  Ri 

c 


i/e 


(32) 


For  the  Vager-Zilitinkavich  model  for  which  Ric  * 1.0 


hr  2.7  — , 
VfN 


(33) 


and  for  the  Mellor-Yamada  model  for  which  Ri  * 0.25 

c 

u 

h - 1.9  — . (3*0 

VfN 

A notable  feature  of  Eq.(32)  is  that  the  mixed  layer  depth  does 
not  depend  upon  the  constant  a which  governs  the  magnitude  of  the  eddy 
coefficients.  Even  for  relatively  small  values  of  0^,  the  mixed  layer 
depth  is  only  weakly  dependent  upon  a.  For  values  of  0 greater  than 
10,  the  slope  of  h*(0^)  is  never  less  than  - 1/6,  and  for  h#  propor- 
tional to  0 "i/9  we  have 
l 
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(35) 


psopav* 


, 1/8  “T 
h * (aRic}  77^71  • 


The  insensitivity  of  the  mixed  layer  depth  to  the  magnitude  of  the  eddy 
coefficients  was  noted  by  Mellor  and  Durbin  (1975)  for  their  integrations 
of  che  Mellor-Yamada  model. 

Hence,  for  wind  dominated  deepening,  the  mixed  layer  depth  is  not 
very  sensitive  to  either  the  magnitude  of  the  eddy  coefficients  or  the 
shape'  of  the  stratification  function  (28).  This  is  evidenced  by  the 
fact  that  the  difference  in  the  critical  Richardson  number  accounts  for 
most  of  the  difference  in  the  mixed  layer  depths  predicted  by  the  Vager- 
Zilitinkevich  and  Mellor-Yamada  models. 

From  an  integral  analysis  of  the  problem  of  wind  dominated  deepen- 
ing of  the  mixed  layer,  Pollard,  Rhines,  and  Thompson  (19,73)  obtained 


h - 1.7  — , 

JfB 


(36) 


which  agrees  fairly  well  with  Eq.  (3M«  the  mixed  layer  depth  predicted 
by  the  Mellor-Yamada  model  for  values  of  the  parameter  0^  in  the  range 
102  - 10s.  In  deriving  Eq.  (36)  Pollard,  Rhines,  and  Thompson  assumed  a 
bulk  Richardson  number  Ri^  for  the  mixed  layer  of  one: 


Pgh(T  - T ) 

RiB  “ — ; — — * 1 » 
® u2  + v2 


(37) 


where  T is  the  mean  temperature  and  u and  v are  the  mean  horizontal 
velocity  components  of  the  mixed  layer  and  T_h  is  the  temperature 
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immediately  below  the  mixed  layer.  It  is  interesting  that  the  integral 
analysis  of  Pollard,  Rhines,  and  Thompson  assuming  a bulk  Richardson 
number  of  1.0  should  give  a result  similar  to  that  of  the  Mellor-Yamada 
model  where  the  criterion  for  the  deepening  of  the  mixed  layer  is  a local 
Richardson  number  of  0.23. 


B.  Heat  Dominated  Mixing 

For  the  case  of  heat  dominated  deepening  of  the  mixed  layer 
with  a constant  wind  stress,  a constant  positive  surface  heat  flux  Q,  and 
negligible  initial  stratification,  the  nondimens ionalization  is  similar 
to  that  for  the  previous  problem  except  that  the  temperature  must  be 
scaled  as  Q/u^  at  . After  the  problem  has  been  scaled,  the  only  variable 
parameter  remaining  is 

0 - . (38) 

2 fu2Ri 

T C 


Hence,  the  fairly  steady  mixed  layer  depth  h achieved  after  several  days 
of  integration  with  a constant  wind  stress  and  surface  heat  flux  is  given 
by 


l 


h 


h*(0  ) • 
2 


(39) 


Figure  9 shows  the  dependence  of  h*  upon  0^  obtained  from  numerical  inte- 
gration of  the  Vager-Zllitinkevlch  and  Mellor-Yamada  models. 

A most  immediately  noticeable  feature  of  Fig.  9 is  that,  although 
the  graphs  of  h*(02)  for  the  Vager-Zilitinkevich  and  Mellor-Yamada  models 
have  a similar  shape,  they  do  not  coincide.  The  difference  in  mixed  layer 


1 

: 

i 


1 
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depth  predicted  by  the  two  models  is  greater  than  that  which  can  bt 

accounted  for  by  the  difference  in  the  value  of  Ri  for  the  two  models. 

c 

Experimentation  with  the  two  models  indicated  that  the  relative  dis- 
placement of  the  two  curves  in  Fig.  9 is  due  principally  to  the  difference 
in  the  value  of  v in  Eq.  (27).  Hence,  the  'shape'  of  the  stratification 
function  (28)  apparently  has  a greater  effect  on  the  mixed  layer  depth 
for  cases  of  heat  dominated  mixing  than  for  cases  of  wind  dominated 
mixing . 

For  increasing  values  of  0 Mellor  and  Durbin  ( 197*0  noted  that 

2 

the  slope  of  h*(0  ) appeared  to  tend  asymptotically  to  - 1.  The  plots 

2 

of  h*(0  ) in  Fig.  9 seem  to  show  the  same  feature.  If  the  dimensionless 
2 

mixed  layer  depth  h*  is  taken  to  be  proportional  to  0-1,  then 


h = 


Ri  uJ 

C T_ 

'ja~  PgQ 


(*<>) 


A notable  feature  of  Eq.  (4-0),  pointed  out  by  Mellor  and  Durbin,  is  that 

the  mixed  layer  depth  does  not  depend  upon  the  Coriolis  parameter  f. 

Hence,  as  the  equator  is  approached  and  the  Ekman  layer  scale  length  u^/f 

tends  to  infinity,  the  mixed  layer  depth  for  heat  dominated  deepening 

becomes  independent  of  f and  dependent  only  upon  the  Monin- Obukov  length 

scale  u3/PgQ. 

T 

For  decreasing  0 the  slope  of  h*( 0 ) increases  and  the  mixed  layer 

2 2 

depth  becomes  more  dependent  upon  the  Ekman  length  scale  and  the  Coriolis 
parameter  and  less  dependent  upon  the  Monin-Obukov  length  scale  and  the 
surface  heat  flux  Q.  For  example,  for  values  of  0^  in  the  range  50  to  100 


40 


k 


(which  corresponds  Co  a range  of  che  surface  heat  flux  of  200  Co  ^00 

cal /an2 -day  when  che  wind  sCress  is  abouc  one  dyne  and  Che  Coriolis 

paramecer  is  abouc  5 x 10  5 sec  x)  che  slope  of  h*(<0  ) is  abouc  - 1/2  and 

2 


h « Ri 


Equation  (40)  is  essencially  the  same  as  the  mixed  layer  depth 
dependence  predicted  by  che  Kraus  and  Turner  model  ( 19^7)  for  heaC  domi- 
nated mixed  layer  deepening  at  all  latitudes.  The  Kraus  and  Turner  model 
assumes  chat  turbulence  within  the  mixed  layer  is  generated  by  wave  action 
rather  than  by  the  mean  velocity  shear  within  the  surface  Etanan  layer, 
and  therefore,  the  depth  of  the  mixed  layer  is  not  influenced  by  the 
earth's  rotation.  Since  the  Kraus  and  Turner  model  assumes  that  th«  rate 

e 

at  which  kinetic  energy  is  imparted  to  the  ocean  by  the  wind  (scaled  by 
Kraus  and  Turner  as  u3)  is  proportional  to  the  rate  of  increase  of  the 

T 

potential  energy  of  the  upper  ocean,  which  is  scaled  as  3gQh  when  the 
mixed  layer  depth  is  fairly  steady  with  time,  we  have 


u3  « 3gQh 
T 


Hence,  the  mixed  layer  depth  dependence  of  the  Kraus  and  Turner  model  for 
heat  dominated  mixing  is  given  by 


Therefore,  a significant  difference  between  the  Kraus  and  Turner  mixed 
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layer  model  and  the  models  of  Vager  and  Zilitlnkevich  and  Mel lor  and 
Yamada  is  that  the  mixed  layer  depths  predicted  by  the  latter  are  inde- 
pendent of  the  Coriolis  parameter  only  in  the  limit  as  the  ratio  of  the 
Elman  layer  length  scale  and  the  Monin-Obukov  length  scale  tends  to 
infinity. 

We  note  here  that  for  an  Ekman  boundary  layer  the  rate  at  which 
kinetic  energy  is  being  imparted  to  the  ocean  by  the  wind  stress  could 
be  scaled  as  the  product  of  the  wind  stress  and  the  mean  velocity  of  the 
Elman  layer*  in  other  words,  as  u*/fh  rather  than  as  u®.  Then,  assuming 
the  rate  of  addition  of  kinetic  energy  to  the  ocean  to  be  proportional 
to  the  rate  of  change  of  the  potential  energy  of  the  upper  ocean,  we 
obtain 


\/fPgQ 


However,  such  arguments  appear  to  give  valid  results  only  for  a small 
range  of  values  of  the  parameter  0^.  Probably  no  simple  theory  of  heat 
dominated  deepening  will  yield  results  that  agree  with  the  results  of 
integration  of  the  Ekman  layer  Eqs.  (22)  over  a wide  range  of  values  of 
the  problem  parameters. 


IX.  SUMMARY  AND  CONCLUSIONS 

The  general  characteristics  of  the  mixed  layer  predicted  by  the 
three  diffusion  models  that  were  compared  are  similar.  Turbulence  gener- 
ated throughout  the  mixed  layer  by  the  instability  of  the  mean  velocity 
shear  keeps  the  layer  mixed  to  a fairly  uniform  density  and  simultaneously 
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creates  a sharp  density  gradient  at  the  base  of  the  mixed  layer  which 


suppresses  turbulent  mixing  and  inhibits  the  diffusion  of  heat  and 


momentum  into  the  region  below 


The  structure  of  the  mean  velocity  field  within  the  mixed  layer 


consists  primarily  of  an  Ekman  spiral  which  is  modified  by  the  stratifi- 


spiral  are  inertial  oscillations  generated  by  recent  changes  in  the 


direction  and  strength  of  the  wind.  Because  of  the  fairly  high  vertical 


diffusivity  within  the  mixed  layer  and  the  large  drop  in  the  vertical 
diffusivity  at  the  base  of  the  mixed  layer,  inertial  oscillations  within 
the  mixed  layer  tend  to  have  fairly  constant  amplitude  and  phase  and  the 
whole  mixed  layer  oscillates  on  top  of  the  thermocline  like  a slab  over 


fication  of  the  upper  ocean,  the  three  diffusion  models  predicted  that 
the  mixed  layer  would  initially  deepen  very  rapidly,  but  that  the  deepen- 
ing would  be  slowed  by  the  Increasing  stratification  at  the  base  of  the 
mixed  layer  after  about  one  inertial  period.  The  Vager-Zilitlnkevlch  and 
Mellor-Yamada  models  predicted  that  the  deepening  of  the  mixed  layer  would 
essentually  stop  after  a couple  of  days,  whereas  the  Munk-A.nderson  model 
predicted  continued  slow  deepening  of  the  mixed  layer  for  the  next  several 
days.  The  difference  in  behavior  of  the  models  is  probably  due  to  the 
fact  that,  for  the  Vager-Zllitlnkevich  and  Mellor-Yamada  models  the  gener- 
ation of  turbulence  by  the  mean  shear  is  completely  suppressed  when  the 
Richardson  number  exceeds  a certain  critical  value  Ri  .whereas,  for  the 


Munk-Anderson  model  the  turbulent  mixing  Is  never  completely  shut  off  as 
the  Richardson  number  increases  and  therefore,  the  turbulence  is  able  to 
erode  the  thermocline  at  the  base  of  the  mixed  layer  even  though  the 
Richardson  number  in  this  region  is  fairly  high. 

A similarity  scaling  of  the  problem  of  wind  dominated  deepening  of 
the  mixed  layer  showed  that  almost  all  of  the  difference  in  mixed  layer 
depth  predicted  by  the  Vager-Zilitinkevich  and  Mellor-Yamada  models  could 
be  accounted  for  by  the  difference  in  the  critical  Richardson  number  of 
the  two  models.  For  values  of  the  dimensionless  parameter  N2Q!a/f®Ri 
corresponding  to  a thermal  stratification  exceeding  10~3  °C/cm  and  a 

-5  -1 

value  of  the  Coriolis  parameter  f of  5 X 10  sec  , the  dependence  of 
the  mixed  layer  depth  h upon  the  problem  parameters  for  the  Vager- 
Zilitinkevich  and  Mellor-Yamada  models  was  found  to  be 

h -2.7  Ri  1/4 , (52) 

C VfN 

where  N is  the  Brunt-Vaisala  frequency  corresponding  to  the  initial 
stratification  of  the  upper  ocean,  u is  the  friction  velocity  associated 

T 

with  the  wind  stress,  and  a is  a parameter  governing  the  magnitude  of  the 

eddy  coefficients.  Since  the  values  of  the  critical  Richardson  number 

Ri  for  the  Vager-Zilitinkevich  and  Mellor-Yamada  models  are  1.0  and 
c 

0.23  respectively,  the  Vager-Zilitinkevich  model  predicted  mixed  layers 
about  30^  deeper  than  those  predicted  by  the  Mellor-Yamada  model  for  the 
cases  of  wind  dominated  deepening  studied. 

For  wind  dominated  deepening,  the  mixed  layer  depth  predicted  by 
all  three  diffusion  models  was  found  to  be  somewhat  insensitive  to  the 
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magnitude  of  the  eddy  coefficients.  Hence,  the  parameterization  of  the 
turbulence  length  scale  for  the  Vager-Zilitinkevich  and  Mellor-Yamada 
models  or  the  determination  of  Kq,  the  vertical  diffusivity  for  unstrati- 
fied turbulent  mixing  for  the  Munk-Anderson  model,  is  not  as  critical  to 
the  prediction  of  the  mixed  layer  depth  as  the  parameterization  of  the 
stabilization  of  the  turbulent  mixing  by  the  stratification. 

For  the  case  of  heat  dominated  deepening  of  the  mixed  layer  with 
a constant  wind  stress,  a constant  positive  heat  flux  at  the  surface,  and 
negligible  initial  stratification  of  the  upper  ocean,  the  three  diffusion 
models  predicted  similar  behavior  for  the  mixed  layer.  The  mixed  layer 
depth  Increased  rapidly  at  first  due  to  the  absence  of  stratification, 
but  then  decreased  suddenly  after  a couple  of  days  as  the  stratification 
in  the  upper  ocean  Increased,  and  stabilized  at  a depth  such  that  the 
stratification  at  the  base  of  the  mixed  layer  was  just  sufficient  to 
stabilize  the  mean  velocity  shear.  The  mixed  layer  then  continued  to  warm 
due  to  the  positive  surface  heat  flux,  whereas  the  region  below  the  mixed 
layer,  where  the  vertical  diffusivity  is  relatively  small,  remained  at  a 
fairly  constant  temperature. 

A similarity  scaling  of  the  problem  of  heat  dominated  mixing  tndi- 

cated  that  the  difference  in  mixed  layer  depth  predicted  by  the  Vager- 

Zilitinkevich  and  Mellor-Yamada  models  was  greater  than  that  which  could 

be  accounted  for  by  the  difference  in  the  critical  Richardson  number  of 

the  two  models.  If  for  the  Vager-Zilitinkevich  and  Mellor-Yamada  models 

the  Richardson  number  is  normalized  by  Ri^,  then  the  Mellor-Yamada  model 

is  found  to  predict  less  suppression  of  the  turbulent  mixing  for  a given 

value  of  Ri/Ri  than  the  Vager-Zilitinkevich  model.  This  difference  in 
c 
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Che  form  of  Che  funcclon  chat  governs  Che  suppression  of  Che  curbulenc 
mixing  by  Che  scracificacion  was  found  co  accounc  for  Che  addicional 
difference  in  mixed  layer  depCh  predicCed  by  Che  Vager-ZiliCinkevich  and 
Mellor-Yamada  models  for  cases  of  heaC  dominaCed  mixing.  Hence,  alchough 
ic  was  ChoughC  chac  perhaps  such  models  could  be  parameterized  jusc  in 
cerms  of  Che  cricical  Richardson  number  characCerisCic  of  che  parCicular 
model,  ic  appears  chac  chis  cannoC  be  done. 

The  similaricy  scaling  of  Che  problem  of  heaC  dominaCed  mixing 
also  indicaCed  chac,  as  noCed  by  Mellor  and  Durbin  ( 197*0 » as  Che  racio 
of  Che  Ekman  layer  lengch  scale  u^/f  and  Che  Monin-Obukov  lengCh  scale 
u^/0gQ  tends  Co  infinity,  che  mixed  layer  depth  becomes  independent  of 
Che  Ekman  layer  lengch  scale  and  dependent  only  upon  Che  Monin-Obukov 
lengch.  Hence,  very  near  Che  equator,  where  che  Coriolis  parameter  f 
tends  to  zero,  Che  mixed  layer  depfh  becomes  independent  of  f. 

As  che  racio  of  che  Ekman  and  Monin-Obukov  lengch  scales  decreases, 

the  Vager-Zilitinkevich  and  Mellor-Yamada  models  predict  Chat  Che  mixed 

layer  depth  becomes  increasingly  dependent  upon  Che  Coriolis  parameter. 

For  example,  for  values  of  0 corresponding  Co  a surface  heat  flux  of 

2 

200  Co  4 00  cal/cm2-day,  a wind  stress  of  about  one  dyne,  and  a value  of 
Che  Coriolis  parameter  of  about  5 X 10  sec  , the  dependence  of  Che 
mixed  layer  depth  upon  Che  problem  parameters  is  approximately  given  by 

1/  u2 

h«Ri  1/4 L.  . (41) 

C ,/fPgQ 

As  Che  surface  heat  flux  Q tends  Co  zero,  it  appears  Chat  Che 
mixed  layer  depth  depends  . only  upon  Che  Elonan  lengch  scale 
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u^/f.  This  is  because,  as  an  initial  condition  for  the  problem  of  heat 
dominated  deepening,  the  upper  ocean  was  assumed  to  be  unstratified. 

Since  the  ocean  is  usually  not  unstratified  to  great  depths,  it  might  be 
expected  that  for  small  values  of  the  surface  heat  flux,  the  depth  of  the 
mixed  layer  would  be  limited  by  the  existing  stratification. 

The  studies  of  the  behavior  of  the  diffusion  models  for  cases  of 
wind  dominated  and  heat  dominated  deepening  of  the  mixed  layer  were  made 
assuming  simple,  idealized  initial  and  boundary  conditions.  In  simula- 
tions of  the  actual  ocean,  the  stratification  is  rarely  ever  linear,  and 
the  surface  wind  stress  and  heat  flux  usually  do  not  remain  constant  for 
extended  periods  of  time.  However,  the  behavior  of  the  mixed  layer  for 
more  realistic  conditions  can  be  generalized  from  the  specific  situations 
studied  here. 

Mixing  within  the  mixed  layer  can  usually  be  classified  as  heat 
dominated,  wind  dominated,  or  convective  depending  upon  the  surface 
boundary  conditions  and  the  stratification.  Given  a specific  wind  stress 
and  nonnegative  surface  heat  flux,  the  mixed  layer  depth  will  tend  towards 
a fairly  stationary  value  where  it  is  in  quasi-equilibrium  with  the  ex- 
isting stratification  and  boundary  conditions.  This  depth  is  limited 
either  by  the  stratification  in  the  case  of  wind  dominated  mixing  or  by 
the  surface  heat  flux  in  the  case  of  heat  dominated  mixing.  If  on  the 
other  hand,  the  net  surface  heat  flux  is  negative,  convective  mixing  will 
be  expected  to  deepen  the  mixed  layer  at  least  to  a depth  sufficient  to 
relieve  the  density  instability  and  probably  somewhat  deeper  due  to  con- 
vective penetration. 
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The  time  required  for  the  mixed  layer  depth  to  adjust  to  a change 
in  the  surface  wind  stress  or  heat  flux  can  be  a few  hours  or  as  long  as 
several  days.  Hence,  if  the  time  scale  of  changes  in  the  surface  boundary 
conditions  is  fairly  short,  as  in  the  case  of  a diurnal  heating  cycle, 
the  mixed  layer  depth  may  seldom  be  in  equilibrium  with  the  surface 
boundary  conditions  and  the  stratification.  As  a result,  the  adjustments 
of  the  mixed  layer  depth  to  short  time  scale  fluctuations  in  the  boundary- 
conditions  may  be  considerably  damped. 

On  a seasonal  time  scale,  the  adjustment  period  of  the  mixed  layer 
is  negligible.  Hence,  the  equilibrium  mixed  layer  depths  determined  for 
the  various  cases  of  mixed  layer  deepening  should  be  useful  for  estimating 
seasonal  mean  values  of  the  depth  of  the  mixed  layer.  For  example,  the 
equilibrium  mixed  layer  depths  determined  for  the  case  of  heat  dominated 
mixing  might  be  used  to  estimate  the  mean  mixed  layer  depth  in  the  spring 
and  summer  when  the  upper  ocean  is  being  warmed  by  solar  radiation. 

There  is  a question  of  how  well  diffusion  models  such  as  those 
compared  here  can  simulate  actual  mixing  in  the  upper  ocean.  This  ques- 
tion is  best  answered  by  comparison  of  the  models  with  observations.  Un- 
fortunately, there  have  been  few  such  comparisons,  mainly  because  of  the 
difficulty  of  obtaining  observations  of  mixed  layer  deepening  suitable 
for  evaluating  the  mixed  layer  models. 

Munk  and  Anderson  ( 19^8)  compared  their  model  with  an  accumulation 
of  data  from  the  Pacific  taken  during  the  summer  at  latitudes  from  15°N 
to  50°N  and  found  that  the  computed  mixed  layer  depths  were  less  than 
half  the  observed  depths.  Munk  and  Anderson  did  not  determine  whether 
the  discrepancy  was  the  fault  of  the  observations  or  the  model.  Denman 
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and  Miyake  (1973)  compared  Che  integral  model  of  Kraus  and  Turner  with 
data  t*ken  at  weather  station  'Papa'  in  the  North  Pacific  and  obtained 
fairly  good  agreement  with  the  observed  temperature  profiles.  Mellor  and 
Durbin  (1975)  used  the  same  data  from  station  'Papa'  to  evaluate  the 
Mellor-Yamada  model  and  also  obtained  fairly  good  agreement  with  Che 
observations.  Since  the  models  of  Kraus  and  Turner  and  Mellor  and  Yamada 
have  some  very  fundamental  differences,  it  might  be  assumed  that  compari- 
sons of  the  mixed  layer  models  with  observations  that  have  been  published 
to  date  have  not  been  sufficient  to  distinguish  one  type  of  model  as  being 
superior  to  another. 

However,  on  the  basis  of  physical  reasoning,  it  does  seem  a severe 
short  coming  of  the  Kraus  and  Turner  and  Denman  models  that  they  do  not 
depend  upon  the  earth's  rotation.  If  the  mixed  layer  is  mixed  by  wave 
generated  turbulence,  then  indeed  the  earth's  rotation  would  probably  not 
influence  the  mixed  layer  depth.  But  there  is  some  doubt  that  mixing  by 
surface  waves  penetrates  sufficiently  to  cause  the  observed  erosion  of 
the  thermocline  at  depths  below  20  to  30  meters.  The  formation  and 
deepening  of  the  mixed  layer  due  to  the  generation  of  turbulence  by  the 
instability  of  the  mean  shear  in  the  wind  driven  surface  Ekman  layer  is 
very  nicely  predicted  by  the  equations  of  fluid  motion  and  simple  tur- 
bulence theory.  And  of  the  arguments  against  this  theory,  that  the  pre- 
dicted mixed  layers  are  not  deep  enough,  that  the  horizontal  length  scale 
is  not  long  enough  for  rotation  to  be  important,  and  that  Ekman  layers 
have  not  been  observed  in  the  ocean,  none  really  holds  up. 

Mellor  and  Durbin  have  shown  that  a stratified,  turbulent,  Ekman 
layer  model  can  predict  mixed  layer  depths  in  close  agreement  with  those 
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that  are  observed.  The  drift  of  the  surface  current  to  the  right  of  the 
wind  stress  as  predicted  by  Ekoan's  theory,  and  the  inertial  oscillations 
of  the  mixed  layer  generated  by  changes  in  the  wind  stress  demonstrate 
the  influence  of  the  earth's  rotation  on  surface  currents  in  the  ocean. 
And  Elonan  spirals  in  the  Great  Lakes  have  been  observed  and  simulated 
using  the  Ekman  layer  equations  by  Gonella  (1971)  and  Csanady  (1972).  It 
should  be  expected  that  the  detection  of  Ekman  spirals  at  sea  might  be 
difficult  due  to  the  problems  of  making  accurate  current  measurements  and 
due  to  the  complicating  effects  of  surface  waves,  inertial  oscillations, 
background  currents,  and  temporal  and  special  variations  in  the  wind 
stress  and  stratification.  Ekman  spirals  may  be  more  frequently  observed 
in  the  ocean  as  current  measuring  techniques  improve. 

Of  the  turbulent  diffusion  models  discussed  and  compared  here,  the 
Mellor-Yamada  model  would  seem  to  be  the  most  suitable  for  modeling  the 
mixed  layer.  The  criterion  of  a Richardson  number  of  about  1 /k  for  the 
stabilization  of  shear  generated  turbulence  by  a density  stratification 
is  a fairly  general  result,  and  it  has  been  seen  that  the  value  of  the 
Richardson  number  for  which  turbulent  mixing  is  suppressed  regulates  to 
a large  extent  the  mixed  layer  depths  predicted  hy  these  types  of  models. 
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